arXiv: 1508.02763v2 [astro-ph.EP] 4 Jan 2016 


Draft version January 5, 2016 

Preprint typeset using style emulateapj v. 01/23/15 


ORBITAL DECAY OF HOT JUPITERS DUE TO NONLINEAR TIDAL DISSIPATION 

WITHIN SOLAR-TYPE HOSTS 

Reed Essick and Nevin N. Weinberg 

Department of Physics, and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 

Cambridge, MA 02139, USA 
Draft version January 5, 2016 

ABSTRACT 

We study the orbital evolution of hot Jupiters due to the excitation and damping of tidally driven 
g-modes within solar-type host stars. Linearly resonant g-modes (the dynamical tide) are driven to 
such large amplitudes in the stellar core that they excite a sea of other g-modes through weakly 
nonlinear interactions. By solving the dynamics of large networks of nonlinearly coupled modes, we 
show that the nonlinear dissipation rate of the dynamical tide is several orders of magnitude larger 
than the linear dissipation rate. We find stellar tidal quality factors Q'^ ~ 10® — 10® for systems with 
planet mass Mp > 0.5Mj and orbital period P < 2 days, which implies that such systems decay on 
timescales that are small compared to the main-sequence lifetime of their solar-type hosts. According 
to our results, there are ~ 10 currently known exoplanetary systems, including WASP-19b and HAT- 
P-36-b, with orbital decay timescales shorter than a Gyr. Rapid, tidally induced orbital decay may 
explain the observed paucity of planets with Mp > Mj and P <2 days around solar-type hosts and 
could generate detectable transit-timing variations in the near future. 


1. INTRODUCTION 

The tide raised by a hot Jupiter excites large amplitude 
waves within its host star. These waves transfer energy 
and angular momentum from the orbit to the star and 
as a result the planet gradually spirals inward. The rate 
of orbital decay is determined by the efhciency of tidal 
dissipation and depends on the amplitude of the waves 
as well as the effectiveness of frictional processes within 
the star. 

Tidal dissipation is often parameterized by the stellar 
tidal quality factor Q'^, where larger Ql implies less dis¬ 
sipation. Perhaps the best constraints on Q'^ for solar- 
type stars come from the observed circulari zation rate 
of solar-type b inaries, which yield Q* 10® (Meibom & 
Mathieu|2005 ). However, because Ql is not a fundamen¬ 
tal property of the star (it depends on the shape and size 
of the orbit and the mass of the perturber), this result 
does not necessarily imply Qi 10® for hot Jupiter sys¬ 
tems. There have been a number of efforts to measure 
Q* from statisti cal modeling of the observ ed sample of 
hot Ju piters (see Ogilvie|2014 for a review). Penev et al. 
(2012) find that the distribution favors Q'^ > 10*^ for a 
specitic set of a ssump tions about the initial conditions. 
Jackson et al. (2008) find a best fit at ~ 10® ® al- 
though they do not rule out much larger values and note 
that it is difficult to obtain tight constraints because of 
the limited sample size and uncertainties in the initial pe¬ 
riod distribution and stellar age. Although there are no 
direct observational measurements of Q* from individual 
hot Jup iter systems (e . g., fro m the detection of orbital 


decay), Jackson et al. (2009) argue that the distribu¬ 
tion shows evidence tor ongoing removal and de struction 
by tides. In addition, Teitler & Konigl (2014) propose 
that the obser ved dearth of close-in p lanets around fast- 
rotating stars (McQuillan et al.||2013) can be attributed 
to tidal ingestion of giant planets. 

Linear tidal driving by the planet resonantly excites 
short wavelength waves within the host star. In solar- 


type stars, these “primary” waves are excited near the 
radiative-convective interface since in this region their 
wavelengths become large and they can couple to the 
long length scale tidal potential. Although the primary 
waves initially have relatively small amplitudes and are 
thus well-described by linear theory, as they propagate 
towards the stellar center their amplitudes increase due 
to geometric focusing (i.e., in order to conserve WKB flux 
within an ever decreasing volume). In hot Jupiter sys¬ 
tems, the primary waves reach large amplitudes as they 
approach the stellar core and become nonlinear, exciting 
many seco ndary waves through nonlinear wave-wave in- 
teractions (Barker fc Ogilvi^|2010 2011 Weinberg et al. 


2012 hereafter WAQB). I'hese secondary waves can have 
much shorter wavelengths than the primary waves and, 
as a result, they can have much larger damping rates 
(due to radiative diffusion). Systems in which nonlinear 
interactions are important may therefore dissipate tidal 
energy much more rapidly than the linear theory esti¬ 
mates. Indeed, in the case of solar-type binaries, the lin¬ 
ear theory estimates yield dissipation rates th at are too 
small by a f actor of > 100 (Q(. ^ 10® — 10^®; Terquem 


et al.||T998| |Goodman fc Dickson||1998| |Ugiivie fc Lm 


2007|). This may indicate that nonlinear processes are 


playing an important role in these systems. 

For a planet with mass Mp > 3Mj(P/day)“®'^ or¬ 
biting a solar-type star, the primary waves reach such 
large amplitudes near the stellar center that the y over¬ 
turn the background stratifica tion and break (Barker| 
'& Ogilvie 2010 Barker 2011). In this strongly non¬ 


linear regime, the primary waves deposit nearly all of 
their energy and angular momentum in a single group 
travel time through the star. The tidal dissipation rate 
therefore equals the energy flux of the initial, linearly 
driven primary waves. The three- dimensional nu meri- 
cal simulations of wave breaking by Barker (2011) yield 


Qi ~ 10®(P/1 day)2-8 for Mp > 3Mj and a solar-type 
star. This corresponds to an inspiral time of « 1 Gyr for 
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a 3Mj planet in a 2 day orbit. 

For a planet with mass 0.5 < Mp/Mj < 3, the primary 
waves do not, in general, break. Nonetheless, they are 
sufficiently nonlinear that they excite many secondary 
waves near the stellar center. In this weakly nonlinear 
regime, the primary waves only deposit a fraction of their 
energy and angular momentum in a single group travel 
time. The value of that fraction, which determines the 
rate of tidal dissipation, depends on the detailed inter¬ 
action between the primary waves and the sea of sec¬ 
ondary waves. The aim of our study is to calculate this 
interaction (and its saturation) in the weakly nonlinear 
regime. Similar types of analyses have been carried out 
in the cont ext of the r-mode instability in spinning n eu- 
tron stars ( Brink et al.||2005 Bondarescu ^al)]|2009 1. 

This paper is structured as follows. In (j we describe 
the formalism we use to study the weakly nonlinear tidal 
interactions and present the equations of motion for our 
mode decomposition. In § |^ we describe how we con¬ 
struct our networks of interacting modes and our method 
for integrating the coupled equations of motion. In §[^we 
present a pedagogical discussion of how different mode 
networks behave. The main results of our calculations 
are presented in § [^ with particular emphasis on the 
tidal evolution of known exoplanetary systems. Finally, 
in §[^we summarize our results and describe some of the 
limitations of our analysis that can serve as directions for 
future work. 


2. FORMALISM 

We are interested in calculating the orbital evolution of 
hot Jupiters due to tidal dissipation within the host star. 
We assume that the planet’s orbit is cir cular, as is the 
case for most of the observed hot Jupiters (Udry & Santos 
2007| |Ogilvie||2014|). If the system is also sufficiently old 
so tha^he planet's rotation is synchronous w ith the orbit 
( Storch fc Lai|20i4 Barker &: Lithwick|2014 ), then there 
is no tidal dissipation within the planet. 

The tide raised by the planet excites a variety of oscilla¬ 
tion modes within the star. Here we limit our analysis to 
solar-type hosts and focus on the excitation of resonant 
g-modes due to linear and (weakly) nonlinear forces. Be¬ 
cause the orbital period of a hot Jupiter is much shorter 
than the rotational period of a solar-type star, the g- 
modes are not strongly modified by Coriolis forces and 
we therefore neglect the star’s rotation. 


is the tidal acceleration, and U is the tidal potential. We 
include only the dominant I = 2 tidal harmonic and since 
we assume that the orbit is circular, 

2 

U{r,t) = -eu;y ^ W2mY2-m{0, (3) 

m —-2 


where e = {Mp/M){R/af, ujq = {GM/Ry/^ is the 
dynamical frequency of a star with mass M and ra¬ 
dius R, Mp is the planet mass, a and H are the orbital 
semi-major axis and frequency, and 14^20 = —(tt/S)^/^, 
W^ 2±2 = (37r/10)^/^, W2±i = 0. We solve Equation Q 
using the method of weighted residuals in which we ex¬ 
pand the six-dimensional phase space vector as 






yr) 


( 4 ) 


where a labels a linear eigenmode with eigenfunction 
eigenfrequency uia, and amplitude qa{t). The sum over 
a runs over all mode quantum numbers and frequency 
signs to allow both a mode and its complex conjugate. 
We normalize the eigenmodes as 

GM^ r 

Eo^ — =2y j (5) 

so that a mode with dimensionless amplitude \qa\ = 1 
has energy Eq- Plugging Equation Q into Equation Q, 
adding a linear damping term, ana using the orthogo¬ 
nality of the eigenmodes leads to a coupled, nonlinear 
amplitude equation for each mode 


(ja + {iiOa +^a)qa = 


Ua{t) -b ^ U*i^{t)qp -b ^ K*^pyq* 


Pi 


( 6 ) 


where 

uy)=-^j d^xpil-VU, (7a) 

uyt) =-^J d^^pi^ ■ {y ■ V) VC/, (7b) 
«a/37 j d^X^^ ■ /2 [1/3, € 7 ] ’ (’^c) 


2.1. Equations of Motion 

We calculate the orbital evolution using the formal¬ 
ism developed in WAQB for studying tides in close bi¬ 
nary systems in which w eakly nonlinear wave interactions 


are im portant (see also Schenk et al. 2001 Van Hoolst 


19941. We now briefly summarize the method and refer 


the reader to WAQB for a more detailed discussion. 

The equation of motion for the Lagrangian displace¬ 
ment ^{r,t) of the stellar fluid at position r and time t 
relative to the unperturbed background is 


= fi K] + /2 [€, €] + pa-tide, (1) 

where p is the background density, fi and /2 are the 
linear and leading-order nonlinear restoring forces, 

atide = - VC/- (I • V) VC/ (2) 


The coefficient 70, is the linear damping rate of the mode, 
Ua and Uaj 3 represent the linear and nonlinear tidal force, 
and Kapj represents the three-mode coupling. 


2.2. Expressions for the Coefficients 

We consider the dynamics of high-order, adiabatic g- 
modes within a solar-type main sequence star. These 
modes are restored by buoyancy and propagate between 
inner and outer turning points determined by the lo¬ 
cations at which uja — N(r) , where N is the Brunt- 
Vaisala buoyancy frequency (|Aerts et al.||2010|). The 
inner turning point is very close to the stellar center 
(?’ct,inner//? — 10“^ (Pct/day)”^) and the outer turning 
point is near the radiative-convective interface at ~ 0.7R. 
Individual modes are described by the quantum num¬ 
bers {I, m, n), where I is the spherical degree, m is the 
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azimuthal order, and n is the radial order. Since the 
modes we consider are all very high-order (n > 50), 
their properties are well approximated by the asymp¬ 
totic WKB expressions given in WAQB. Using a 5 Gyr 
old solar model from the EZ code (Paxton]2004 ), we find 


UJa CS 7-Wo, 

n„ 


7„ cs 2 X 10-“A2 



(8a) 

(8b) 


where A^ = la{la + !)• The dominant linear damping 
mechanism of the high-order g-modes is radiative dif¬ 
fusion of the temperature fluctua tions that accompany 


the mode density perturba tions (Terquem et al. 1998 


Goodman & Dickson 19981. Shorter wavelength modes 


therefore have larger damping rates. 

By plugging Equation ra into Equation (7a), we can 
express the linear driving coefficient Ua in terms of the 
dimensionless linear overlap integral 


/a = ^ / d^^pC ■ V (9a) 

~ 2.5 X 10-3 , (9b) 


where the numerical expression assumes la = ‘2, rua = m 
and is accurate for mode periods Pa > 0.3 day (Figure 
11 of WAQB). Low-order, I = 2 g-modes have large la 
but UJa ^ they comprise the quasi-static response 
of the fluid (the equilibrium tide). High-order, I = 2 
g-modes have small la but can nonetheless have large 
linear amplitudes if Wq ~ 2D; they comprise the resonant 
response of the fluid (the dynamical tide). 

The three-mode coupling coefficient Kapj is symmetric 
under the interchange of mode indices. Angular momen¬ 
tum conservation leads to the following angular selection 
rules for the three modes: (i) la + 1/3 + Ij must be even, 
(ii) nia + mp + = 0, and (hi) the triangle inequality, 

\la — lp\ ^ l'^ la+l/s- We focus On the parametric insta¬ 
bility involving three-mode interactions between a high- 
order “parent” g-mode and a pair of high-order “daugh¬ 
ter” 5 -modes whose summed frequency nearly equals the 
parent’s frequency. For such a triplet, the coupling is 
strongest in the stellar core, where the Lagrangian dis¬ 
placements of the modes peak. For a solar-type star (Ap¬ 
pendix A in WAQB) 


tide is small by comparison. We therefore restrict 
our mode networks to parent modes that comprise the 
dynamical tide response of the star (i.e., linearly resonant 
parents) and ignore the equilibrium tide response and 
nonlinear tidal driving. 

Finally, we assume that only linearly resonant modes 
(parents) have non-zero linear tidal forcing Ua- This is 
justified because the linear forcing coefficient Ua is much 
smaller for daughter modes and their driving is far off res¬ 
onance. Its secular effect will therefore be negligible com¬ 
pared to the resonant three-mode interactions. Ignoring 
such forcing allows us to adopt a convenient change of 
coordinates that significantly spe eds up the integration 
of the amplitude equations (see § 3.3). 


2.3. Instantaneous Orbital Decay Time-scale 

The energy of the stellar modes is E^/Eq =J2 Qa^a + 
{l/3)J2ka/3-yiqa<l/3<l-y + c.c). The rate of energy loss 
within a solitary star is therefore 


E*. 

Eq 


E* 


X! (daQa + C-C) + ^ kap^ (90,9/397 + C.c) 

O' a^7 

-2 ^ 7a9a9a + kap-y'Ja (9q 9/397 + C-c) , 

O' OL^-) 

-2^JaEa, ( 11 ) 

O' 


where we substituted the equations of motion (Equation 

S for the mode amplitudes and neglected the terms from 
e three-mode couplings because they are much smaller 
than Ea = 9 q 9 q£'o. This is the rate at which energy is 
dissipated within the star by radiative diffusion. 

The dissipation of tidally excited stellar modes removes 
energy from the orbit, and the orbit therefore decays. We 
assume that the only dissipation in the system is due to 
the linear damping of waves excited within the star. Al¬ 
though the rotational energy of a synchronized planet in¬ 
creases as the orbit decays, this change is small compared 
to the corresponding change in orbital energy. Similarly, 
the energy in the excited stellar modes themselves may 
change with orbital period, but this also is a small effect 
(see Appendix . 

Because |F*/Forb| ^ H, where Forb = —GMMp/2a 
is the orbital energy, we model the back-reaction on the 
orbit as a steady decrease in Forb of quasi-Keplerian cir¬ 
cular orbits. The timescale of the instantaneous, orbital 
energy decay is then given by 




~ 2 X I03 


T 

0^ 


1 day 


( 10 ) 


where Pa is the period of the parent mode and T « 
0.1 — 1 is an angular integral that depends on each mode’s 
I and m. The coupling occurs mostly near the parent’s 
inner turning point Ta^inner and scales as P^ because the 
parent’s displacement there varies as ^a ^g inner ^a- 
Although the equilibrium tide amplitude is large, its 
three-mode coupling cancels signific antly with nonlinear 
tidal driving Uap (WAQB; see also Venumadhav et al. 
2014). As a result, for a hot Jupiter system, the nonlinear 
dynamics are dominated by three-mode coupling to the 
dynamical tide; the energy dissipated in the equilibrium 


Aorb 



( 12 ) 


at each P, and we can compute a corresponding time- 
averaged decay time-scale 


(te) = 


a 


-^orb 

wy 


(13) 


where (E^) is the time-averaged energy dissipation rate 
with the average spanning several resonance peaks (if 


^ Turbulent dissipation of the e nuilibrium tide within t he con¬ 
vection zone yields Ql 10^ ~ 10^ (jPenev Sasselov|2011 1 , which 
is much larger than the we find due to nonlinear daniping of 
the dynamical tide (§[^. 
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{ te ) oc a" then the “inspiral time” into the star will be 
{TE)/n). We describe our method of time-averaging in 
Appendix Using the language of linear tidal theory, 
{te) is often parameterized in terms o f the st ar’s tidal 
quality fact or (|Goldreich & Soter|| 19661 see also | Jackson 
eTallpOSl ) 


Q* = 7.5 X 10' 


,6 [ {jA 

\Gyr 


Mp 

Mj 


P 

day 


- 13/3 


(14) 


where the expression assumes a circular orbit about a 
solar-type star and ignores dissipation within the planet. 
Although is often taken to be a constant and fun¬ 
damental property of the body, in general it depends 
on the companion mass, orbital frequency, and the tidal 
harmonics (l,m). 


3. BUILDING AND INTEGRATING THE MODE 
NETWORKS 

In the absence of nonlinear three-mode interactions, 
the energy of a linearly driven parent mode a is (Equa¬ 
tion 29 in WAQB) 


Alin 

Eq 






(15) 


where Aq, = Wq, —mafl is the linear detuning. The parent 
is unstable to nonlinear three-mode interactions if there 
exists a pair of daughter modes j3 and 7 such that E\\-n > 
Athr, where the threshold energy is (see Appendix [b| 

Athr _ 1 f Ibl-t \ 


1 + 


7 / 3+77 


(16) 


Here = top + -\- nia^ is the nonlinear detuning 

of the daughter pair. Daughter pairs with smaller |A^-y| 
and larger |Ka/ 37 l (he., stronger nonlinear coupling) yield 
smaller Athr and are more readily unstable. In a three¬ 
mode system, unstable daughters with small initial am¬ 
plitude undergo a phase of exponential growth at a rate 


hsind ~ |lT/a/37l ^ 0 - 


(17) 


Eventually, the daughters reach an energy comparable 
to or greater than the paren t’s a nd the system reaches a 
nonlinear equilibrium (see § 4.2 and Appendix . 

Eor the tide raised by even a 0.1 Mj companion in a 
3 day orbit, ther e ar e ^ 10^ daughter pairs for which 
Athr < Alin (see §[^. In § El we systematically explore 
the dynamics of large multi-mode, multi-generation sys¬ 
tems. In brief, we find that the parent drives many of the 
unstable daughters to large amplitudes and these daugh¬ 
ters, in turn, drive granddaughters to large amplitudes, 
and so on. The total number of potentially unstable 
modes and the number of couplings is larger than the 
number we can integrate on a computer in a reasonable 
time (~ 10^ and ^ 10^, respectively). The issue then 
is whether we can reliably calculate the total tidal dis¬ 
sipation rate with a mode network that contains only a 
subset of the potentially unstable modes. We will present 
evidence in § |^ that this is possible but we must build 
our networks carefully and systematically. 

In §§ |3.1| and |3.2[ we describe how we build networks 
consisting of sets of three-mode couplings (i.e., sets of 
triplets) and collective couplings, respectively. And in 



Figure 1. Cumulative distribution of -Fthr for = O.lMj and 
P ~ 3 days (257928 sec), (blue) Couplings from parents to daugh¬ 
ters. The vertical blue line corresponds to the parent’s linear en¬ 
ergy Piin- (red) Couplings from daughters to granddaughters. 


§ 3.3 we describe our method for integrating the coupled 
rnode amplitude equation (Equation^. 

3.1. Building Three-mode Networks 

Although there are many daughter pairs with Athr < 
Alin, we show in § 13 that pairs with low Athr domi¬ 
nate the dynamics olTarge multi-mode systems. We find 
that if we gradually increase the size of our networks by 
adding pairs with progressively higher Athr, the system 
converges to a dissipation rate A* that does not change 
signihcantly as we add even more modes. We must also 
include a sufficient number of generations (at least par¬ 
ents, daughters, and granddaughters) in order to obtain 
convergent results. Therefore, to build our mode net¬ 
works, we comprehensively search the mode parameter 
space and construct, for each generation, a complete list 
of pairs ranked by Athr- 

In order to carry out our s earch, w e us e th e expressions 
for w, 7 , and k (Equations 8 a 8 b and 10) to solve for 
Athr- For a given parent mode a, we first nnd the local 
minima of Athr in the daughter parameter space {{np, Ig, 
m/s), (uj, l-y, m-y)}. In general, Athr is minimized approx¬ 
imately where the sum in quadrature of A^y and 7 ^ -|- 7 y 
is minimized (modulo the angular selection rules and a 
relatively weak dependence on the angular integral A). 
Daughters with higher I have smaller A,gy (because they 
are more densely spaced in frequency) but larger 7 (be¬ 
cause 7 ~ P); the regions of small Athr therefore occur 
where these two countering effects are balanced. After 
finding the local minima, we expand our search around 
those minima and find pairs with progressively higher 
Athr- Because 7 ^ ^t high enough I the damping 
dominates detuning, and Athr increases with increasing 
1. We truncate our search upon reaching an Zmax such 
that Athr > Ahn (i-C., a stable triplet). In practice, we 
find that for parent-daughter coupling, the dissipation 
A* is dominated by the 10-100 lowest Athr triplets. 

Figure shows the distribution of Athr for parent- 
daughter coupling assuming a 0.1 Mj companion in an 
orbit near three days. There are a few pairs with very 
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low i?thr because they just happen to have particularly 
small despite having I < 3. There is a much larger 
sample of ^ 10 ^ pairs that have larger A^-y and/or 7 
which still yield ifthr "C -Biin- 

We carry out a similar search when we consider the 
coupling of daughters to granddaughters. We show the 
distribution of ifthr for daughter-granddaughter coupling 
in Figure The ifthr of the most unstable daughter- 
granddaugE^ter triplets is much smaller than the ifthr of 
the most unstable parent-daughter triplets (i.e., the red 
curve in Figure is far to the left of the blue curve). 
Because Kap-y oc uj~'^ and AjSy/uja oc Ua, for low Ethr 
pairs we find ifthr oc {Apy/KapyOJaY oc w®. The factor 
of two decrease in frequency with each generation there¬ 
fore means that ifthr decreases by ~ 2 ® (for a full discus¬ 
sion, see Appendix]^. Physically, Athr decreases because 
lower frequency mooes (i) penetrate deeper into the core 
where Lagrangian displacements are larger, and (ii) are 
more densely spaced in frequency and therefore can have 
smaller detunings. As a result, each generation is ev¬ 
ermore susceptible than its predecessor to three-mode 
instabilities. This has important implications for the dy¬ 
namics of large multi-mode, multi-generation systems, as 
we describe in § |^ 

3.2. Building Collective Networks 

For high-order g-modes, the frequency spacing between 
neighboring modes is |Au;| ~ oj/n ^ w. Therefore, if a 
pair of daughters is resonantly excited by a parent, there 
is a good chance that neighboring modes will also be 
resonantly excited by that same parent. The dynamics of 
such a system can be very different from that of a simple 
three-mode system; in particular, the daughters can grow 
as a single, collective unit with growth rates that are 
much higher than the three-mode case (see WAQB and 
Appendix . 

To appreciate why collective sets can grow so quickly, 
consider a simplified system in which a single parent 
mode a is coupled to N daughter modes that are closely- 
spaced neighbors in (l,n) space. A study of the dynam¬ 
ics of such a system reveals that the modes all oscillate 
nearly in phase with each other. The equations of motion 
for each of the N daughter modes can thus be approxi¬ 
mated as 


qp + {iujp + ^p)qp = iu}p i^lcPyCdy 

7 

~iu}pNK*^p^q*^q*. (18) 

The dynamics look like the three-mode case, but with an 
effective coupling coefficient that is N times larger. In 
particular, the instability growth rates (threshold ampli¬ 
tudes) are approximately N times larger (smaller) than 
the three-mode case. 

When building our mode networks, we use separate 
algorithms to search for collective sets and three-mode 
sets. A simple but incomplete way to build collective 
sets is to first find a daughter with a frequency nearly 
equal to half that of the parents and then progressively 
add neighbors with An = ±1, ± 2 , ±3,.... At first, ifthr 
will decrease as more modes are added and N increases. 
However, for large enough An, the detuning of the outer 
most modes becomes so large that adding more modes 


does not decrease ifthr any furtherR More detail is pro¬ 
vided in Appendix |E.4[ Although this method naturally 
picks out collective sets (and is similar to the approach 
described in WAQB), it potentially misses many collec¬ 
tively unstable modes. For example, there can be distinct 
groups of modes that are not nearby neighbors and yet 
together form a collective set. For this reason, we use a 
more sophisticated method when building collective net¬ 
works. We describe this method in Appendix [P] 

In § 1^ we show that collective sets are excited and 
initially grow much more rapidly than three-mode sets. 
However, when the entire network ultimately reaches its 
nonlinear equilibrium and saturates, we find that the col¬ 
lective sets do not alter A* significantly. We therefore 
find that we can accurately calculate if* with networks 
that include only three-mode sets. 


3.3. Integration Method 

We integrate the amplitude equation (Equation]^ for 
each mode of a network using an adaptive step-size 4*^- 
5*^ order Runge-Kutta integrator. Our integrations take 
advantage of a convenien t chan ge of coordinates, also 
described in Brink et al. (2005). The integration step 
size is limited by the fastest frequency in the equations. 
Because the linear and nonlinear forcings all involve res¬ 
onant interaction^ the linear and nonlinear detunings 
(Ac and Apy) are all small (<C H). In fact, the fastest 
time scale in these equations is typically the natural fre¬ 
quency bJa of each mode. By changing coordinates to 
Xa = qaC‘^°‘*i we can remove these frequencies from the 
equations of motion at the cost of adding a slowly vary¬ 
ing time-dependent term to each three-mode coupling. 
This increases the typical integration step size by approx¬ 
imately the ratio of Wq to the detuning (« 10 ^ — 10 "^). 

In order to further speed-up the integrations, we par¬ 
allelize across multiple CPUs. We achieve this using 
standard parallelization techniques, with care taken to 
equally distribute the amount of work across each CPU. 
For example, the computation of Xa scales with the num¬ 
ber of couplings included for that mode. Therefore, when 
we parallelize the computation of Xa by splitting modes 
among processes, we attempt to divide modes into sets 
with equal numbers of couplings, rather than equal num¬ 
bers of modes. We test several different parallelization 
methods, including an implementation using Python’s 
subprocess module. Python’s multiprocessing module, 
and a Python wrapper for OpenMPI. All our implemen¬ 
tations scale better than A^cp^, although which imple¬ 
mentation is fastest depends on specifics of the hard¬ 
ware. The Python multiprocessing implementation gen¬ 
erally performed best, and parallelized across 15 2.7 GHz 
Quad-Core AMD Opteron Processors, it takes ~ 40 sec¬ 
onds to integrate one of our largest networks (~ 2.4 x 10"^ 
modes with ~ 3 x 10® couplings) through 10 orbital pe¬ 
riods. 


4. MODE DYNAMICS 


^ In addition, the magnitude of Kapy becomes small for An > Ua 
because the coupled daughters are no longer spatially resonant with 
the parent (see Figure 12 in WAQ B). 

^ For reasons described in § |2.2| we assume that only the linearly 
resonant parent modes have a non-zero linear tidal forcing Ua. 
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Figure 2. Energy dissipation rate Et as a function of orbital period P for networks that include different numbers of mode generations 
Ngena ■ We show results for Mp = Mj and orbits near P = 1 day (left panel) and P = 3 days (right panel). The sharp linear resonance peaks 
occur when the tidal driving frequency is resonant with an / = 2 g-mode of the star, (blue squares) Networks with only the ten most resonant 
parents, (red circles) Networks with parents and daughters, (green triangles) Networks with parents, daughters and granddaughters (shown 
here are the results of our “reference network”; see § |4.6| l . (gray lines) The analytic estimate of the steady state dissipation for parent-only 
networks that contain t he 1 00 most resonant parents, (red crosses) The analytic estimate of the steady state dissipation for parent-daughter 
networks (see Equation |21[|. We do not have an analytic estimate for parent-daughter-granddaughter networks. 


In order to build up intuition for the results from large 
multi-mode, multi-generation networks, we describe the 
mode dynamics of increasingly complicated networks. 
We begin in § |4.1| with a network that consists only of lin¬ 
early resonant parents (i.e., we ignore all nonlinear cou¬ 
plings) and show that our simulations rec over the dissi¬ 
pation rates of standard linear theory. In § |4.2| we couple 
linearly resonant parents to unstable daughter modes but 
do not allow the daughters to couple to granddaughters. 
We find that including even just this first generation of 
nonlinear couplings enhances the dissipation rate by a 
factor of « 100 (« 10) relative to the lin ear result for 
P = 1 day (3 day) and Mp = Mj. In § 4.3 we allow 
the daughters to couple to granddaughters and find that 
this further enhances the dissipation, yielding a rate that 
is « 10® (« 10®) times larger than the linear r esult for 


P = 1 day (3 day) and Mp = Mj. We find in §§ 4.4 and 


|4.5| that the dissipation rates do not change significantly 
when we include even more generations (great grand¬ 
daughters and beyond) and collective sets, respectively, 
suggesting that th e sys tem has reached a convergent, sat¬ 
urated state. In § |4.6| we explore the minimum network 
size needed to attain such a convergent state. 

4.1. Linear Parents Only 

If we include only linearly driven parents in the net¬ 
work, then \E^, \ = 2 ^ 7^ (illc()iin, where the parent linear 
energy (illct)iin is given by Equation (151. The dissipa¬ 
tion is typically dominated by the most linearly resonant 
parent, although other modes can contribute if no single 
mode is particularly resonant. In Figure we show E* 
due to the ten most resonant parents over a small range in 
orbital period. In the absence of nonlinear interactions, 
the orbit evolves rapidly through the sharp resonance 
peaks where is large. As we describe in Appendix 
[A} the time average dissipation rate is the sum of the 
instantaneous E* weighted by the amount of time spent 


at that period (see alsq_ Goodman & Dickson 1 998| ). In 
the left panel of Figure we show te (Equation |T^ due 
to the single most resonant parent assuming P ~ 3 day , 
Mp = Mj. An analytic calculation using Equations (8b I 
and (15), and assuming uJal^Ua ^ 7a, yields (see 

Appendix 


{TE)rm — 1-4 X 10 


which translates to 


12 


Mj) Vday 


yr, 


1,1 X lO-" ( P 


- 4/3 


(19) 


( 20 ) 


This is in good agreement with our numerical integra¬ 
tions. The Mp dependence of {TE)iin is due to the linear 
forcing coefficient [/„ and the dependence on P is due to a 
combination of 7a, Aa, and Ua- We will show that when 
we include nonlinear interactions, the instantaneous de¬ 
cay time has a dramatically different magnitude and scal¬ 
ing with Mp and P. 

The right panel of Figure shows the results for the 
same parameters as the left panel, but for networks with 
either the 10 or 25 most resonant parents. Despite in¬ 
cluding these additional modes, the instantaneous decay 
time is nearly identical in both cases. More generally, 
we find that including multiple parents has very little 
effect on the total dissipation (even for networks with 
nonlinear interactions) as long as P < 4 days. This is 
because, for P < 4 days, the parent mode spacing is suf¬ 
ficiently sparse that the most resonant parent typically 
has a much larger E\in than the neighboring parents, and 
it therefore dominates the dynamics (Appendix [C|). 


4.2. Parents and Daughters 

In Figure 15 we show the mode dynamics of networks 
that include daughters (but not granddaughters) coupled 
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Figure 3. Orbital decay timescale te as a function of the number of modes AlmoUgg for networks that include different numbers of 
mode generations Algens- We show results for a Jupiter-mass companion orbiting at a period near three days (257928 sec, chosen to be 
approximately half way between a resonance trough and resonance peak for our stellar model). The left panel networks have a single 
parent mode and the right panel networks have 10 or 25 parent modes, (blue squares) Networks with only parents, (red circles) Networks 
with parents and daughters, (green triangles) Networks with parents, daughters, and granddaughters. The filled triangle corresponds to a 
reference network with collective granddaughter modes added, (purple pentagons) Networks with parents, daughter, granddaughters, and 
great-granddaughters. The structure of the networks with only a single parent (left panel) are as follows: fVgena = 2 networks range from 
one daugh ter p air up to 500 daughter pairs. Ngena = 3 networks mostly correspond to 10 daughter pairs and either 25, 50 (our reference 
network; § |4.6[ l, 75, 150, or 200 granddaughter pairs per daughter, although there are networks with 50 and 200 daughter pairs, each with 
50 granddaughter pairs per daughter. Algena = 4 networks are all extensions of our reference network, adding 10, 25, 50, 100, 200, or 500 
great-granddaughter pairs per granddaughter. The structure of the multi-parent networks (right panel) are as follows: either 10 or 25 
parent modes; 0, 10, 25, or 50 daughter pairs per parent; and either 0, 50, 100, or 200 granddaughter pairs per daughter. 


to a linearly resonant parent. The top panel shows a 
simple three-mode system involving a parent coupled to 
only its lowest i?thr daughter pair. Initially, the daugh¬ 
ters are at small energy and the parent is at its linear 
energy Eym- Because Enn > it'thr, the system is unsta¬ 
ble and the daughters underg o a rapid initial growth at 


the rate given by Equation (171. Eventually the sys¬ 
tem reaches a nonlinear equilibrium in which the parent 
has energy Ea = Ethr and the daughters have energy 
Ep^j ~ \Ua/2Ka/3’y\Eo (see Appendix [ b] and WAQB). 

The middle panel of Figure shows the same par¬ 
ent now coupled to the ten lowest Athr daughter pairs. 
Because Apn > Ethr for all ten pairs, initially all the 
daughters grow. However, eventually the parent energy 
drops to the minimum Ethr and only the lowest thresh¬ 
old daughter pair remains excited; the other daughters 
decay due to linear damping. The nonlinear equilibrium 
of this system is therefore equivalent to the three-mode 
network shown in the top panel. 

The network shown in the middle panel of Figure [^as¬ 
sumes that each of the ten triplets only share a parent. If 
the triplets also share daughters (e.g., daughter a couples 
to daughter b and daughter c), then the dynamics can be 
more complicated. Such a network is shown in the bot¬ 
tom panel of Figure^and illustrates how such a dditional 
couplings parasit ically excite other daughters O’Leary 
Ifc Burkart]l2ni4l consider a similar mechanism in order 
to explain the odd resonances observed in the KOI-54 
light curve. Despite these additionally excited modes, 

This is a form of nonlinear inhomogeneous driving and is de¬ 
scribed in WAQB. 


the lowest Ethr pair still dominates the dissipation. 

In Figure we show te for networks that include 
only parent-daughter couplings (Agen = 2) assuming 
E ~ 3 day, Mp = Mj. We find that at this period parent- 
daughter coupling decreases te by a factor of ~ 10 rela¬ 
tive to the linear result. Numerically, both te and (te) 
are nearly independent of the number of daughter modes 
in the parent-daughter networks because the daughter 
pair with the lowest Ethr dominates the dynamics and 
dissipation. The other daughters, while excited, do not 
reach significant amplitudes and therefore have little ef¬ 
fect. Even for large numbers of modes, parent-daughter 
systems behave much like those shown in the middle and 
bottom panels of Figured We therefore find that (te) 
is well-approximated by trie analytic calculation that as¬ 
sumes only one parent and its single lowest Ethr daughter 
pair (see Appendix [G| and Figure [To|) 

/ p \ 19/6 

{TE)p.d ^ 2.0 X 10“ f — j yr, (21) 


assuming I = 1 daughters and 




1.5 X 10® 


\Mj 



( 22 ) 


The agreement between the numerical result for parent- 
daughter networks containing multiple daughters and the 
above three-mode estimate is further illustrated in Fig¬ 
ure The open circles show the numerically computed 
E* of a parent-daughter network consisting of the ~ 20 
lowest Ethr daughter modes. Each of these open circles 
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Figure 4. Parent-daughter networks (A^gens = 2) with different 
structures. The left panels show the energy Ea of each mode as 
a function of time. In each of these panels, the parent is the blue 
line and is initially at Eym but ultimately settles into a nonlinear 
equilibrium at an energy E <C Ey^. The right panels show the 
coupling diagrams in the n-l plane, with circles representing the 
included modes and line connections indicating the coupling struc¬ 
ture. (top) A parent and the lowest .Ethr daughter pair, (middle) A 
parent and the ten lowest .Ethr daughter pairs with each daughter 
mode couple to only one other daughter mode, (bottom) A parent 
and the ten lowest -Ethr pairs with all allowed couplings between 
daughters. 

is covered by an “x”, which represent the analytically 
computed assuming only the minimum iiithr daugh¬ 
ter pair. 

We will now see, however, that the dynamics are much 
more complicated when granddaughters are included, 
with many more modes excited to significant amplitudes. 

4.3. Parents, Daughters, and Granddaughters 

In the absence of granddaughter couplings, the lowest 
Ethr daughter pair settles into a nonlinear equilibrium at 
an e nerg y Ep^j ~ \Ua/2Kai3^\Eo. However, as discussed 
in § |3.l] there are many granddaughter pairs that are 
unstable to such high energy daughters (see Figure [^. 
The parent-daughter solutions of the previous section are 
therefore unstable and never realized. 

For very small networks that include granddaughters 
we sometimes observe periodic limit cycles. However, for 
even slightly larger networks with more complicated cou¬ 
pling topologies, the limit cycles begin to take on a more 
chaotic appearance. And for the very large networks 
that we find yield convergent dissipation results (> 10^ 
modes), the dynamics cease to display any clear limit 
cycle behavior over long time scales and instead show 
persistent large amplitude fluctuations involving many 
excited modes (Figure [^. 

We can roughly unaerstand the behavior of these 
networks using intuition from simple three-mode sys¬ 
tems. Initially, an unstable, linearly driven parent excites 
daughters to large energy. The daughters drain energy 
from the parent and the parent’s energy drops. How¬ 
ever, the daughters then excite granddaughters and the 
daughters’ energy drops. The daughters no longer drain 
enough energy from the parent and the parent begins to 
recover due to linear driving. The rising parent excites 



10-20 10-19 10-18 




Figure 5. Cumulative distributions of the number of modes in¬ 
cluded in the network (blue) and the energy dissipation rate of the 
modes (red) as a function of each modes .Ethr- For granddaughter 
modes coupled to more than one daughter, we take their minimum 
Ethr- (top) ^ network consisting of one parent, 50 daughter pairs, 
and 50 granddaughter pairs per daughter, (bottom) a network con¬ 
sisting of one parent, 500 daughter pairs, and 50 granddaughter 
pairs per daughter. Approximately 90% of the energy dissipation 
is due to the first 50% (40%) of modes ordered by Ethr for the 
smaller (larger) network. 


the daughters again and the cycle restarts. 

Unstable granddaughters have lower frequencies and, 
in general, higher I than the parents and daughters. They 
therefore often have much smaller radial wavelengths 
(i.e., much larger n) and thus much larger linear damping 
rates (7 oc n^). This means that granddaughters can dis¬ 
sipate energy more rapidly than daughters even if they 
are at a lower amplitude. 

In Figure]^ we show te for large networks that include 
parents, daughters, and granddaughters (Vgen = 3) as¬ 
suming P ~ 3 day and Mp = Mj. We see that networks 
with granddaughters are more dissipative than parent- 
daughter only networks and yield te < lO^^rejin at 
P ~ 3 day. The figure also shows te as a function of 
the number of modes Amodes in the network. We find a 
systematic uncertainty in te associated with the struc¬ 
ture of the network. However, this uncertainty is small 
compared to the increase in dissipation associated with 
the inclusion of granddaughter modes. In this sense, we 
find that te is not particularly sensitive to the number of 
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granddaughters nor to the details of the network struc¬ 
ture as long as the number of granddaughter modes is 
sufficiently large (> 10^). We illustrate this po int f ur- 
ther when we discuss our reference networks in § |4.6[ 

When we build networks with larger Wnodes) we do s o 
by adding modes of increasingly larger iSthr (see 13.1). 
The fact that te does not change as we increase Wiodes 
above ^ 10^ suggests that the lowest Ethr modes dom¬ 
inate the energy dissipation and, therefore, that our 
method for building networks reliably captures the bulk 
of the dissipation. We illustrate this more explicitly in 
Figure which shows that the overwhelming majority 
of the energy is dissipated by the modes with the low¬ 
est Ethr and that modes with ever larger Ethr contribute 
less and less to the total dissipation. This suggests that 
selecting modes based on their Ethr identifies the dy¬ 
namically relevant couplings and that including enough 
modes in this way yields convergent results. 

4 . 4 . Great Granddaughters and Beyond 

If daughters and granddaughters are excited, what 
about great granddaughters (iVgen = 4) and so on? This 
s eems particularly likely given that Ethr oc (see 

§ 3.1). Indeed, based on our experiments with networks 


that include up to five generations, we observe that the 
cascade continues into many generations. However, as 
long as we include enough modes, we find that the total 
dissipation rate plateaus once we include granddaugh¬ 
ters. In effect, we do not need to resolve the innermost 
scales of the energy cascade in order to obtain an accu¬ 
rate estimate of E^. 

We illustrate this in the left panel of Figure which 
shows Te for networks that go up to A^gen = 4. There are 
dramatic decreases in te when going from just parents 
(fVgen = 1) to parents and daughters (fVgen = 2), and 
again when adding granddaughters (iVgen = 3). However, 
we see only a slight decrease in te when we add great- 
granddaughters (iVgen = 4). In particular, for sufficiently 
large iVgen = 4 networks (iVmodes ^ 2 x 10^), we find that 
Te plateaus at a value that is only ~ 3 times smaller than 
that of our iVgen = 3 reference network at this P. This 
suggests that truncating at iVgen = 3 yields reasonably 
accurate estimates of te- 


4 . 5 . Dynamics of Gollective Networks 

WAQB showed that sets of modes can be collectively 
unstable even if each pair within the set is stable by it¬ 
self. We describe algorithmi c app roaches to identify and 
select such sets of modes in § |3.2| and Appendix [E These 
collective sets can have growth rates that are hundreds of 
times fast er t han the separate three-mode growth rates 
(Equation 17). We illustrate this in the left panel of Fig¬ 
ure whicn shows several collectively unstable sets of 
modes being rapidly driven to large amplitudes after only 
a few hundred orbital periods. However, for P < 10 days 
(see Figure 7 of WAQB), the lowest Ethr for individual 
three-mode triples is lower than the collective stability 
threshold. This means that after the collective modes 
grow rapidly, the parent’s amplitude is still large enough 
to drive three-mode triples. The right panel of Figure 
shows that the slowly growing three-mode triples even¬ 
tually reach large amplitudes and drive the parent below 
the collective instability threshold. At that point, all the 



Figure 6. Parent-daughter networks (A^gens = 2) that include 
collective sets in addition to three-mode sets, (black dashed line) 
the parent mode, (black solid line) daughters selected with three¬ 
mode algorithm, (grey solid lines) daughters selected with collec¬ 
tive algorithms. In the left panel we show the dynamics of the 
system at early times. Several separate sets of collective modes are 
excited, grow rapidly, and reach significant energies by £ ps 200P. 
In the right panel we show the dynamics out to much longer times. 
The slowly growing set of unstable three-mode daughters eventu¬ 
ally reach large amplitudes at £ ps 3x 10‘^P. This drives the parent’s 
amplitude down and the collective modes decay away. We only plot 
one out of every ten collective modes from this simulation. 


collective modes “turn off” and decay, leaving the steady 
state predicted by simple three-mode systems. 

The network in Figure only includes parents and 
daughters. However, because these d augh ters are unsta¬ 
ble to granddaughter interactions (§ 4.3), the collective 
modes may not decay forever but instead may saturate 
at non-trivial amplitudes. In principal, because collective 
modes can have significantly larger I than the minimum 
Ethr pair and thus larger damping coefficients, they may 
dissipate energy more rapidly. Nonetheless, numerical 
experiments reveal that sufficiently large networks con¬ 
structed out of only three-mode pairs yield nearly the 
same A* as networks that also include collective excita¬ 
tions. This can be seen in Figure the fille d triangle 
corresponds to the reference network (§ 4.6) with the 
addition of collective granddaughters. Because collective 
networks are expensive to simulate and do not change 
the calculated A*, from here on we do not include them 
in our calculations. 


4 . 6 . Reference Network Integrations 

To summarize, we find that networks with parents, 
daughters, and granddaughters yield convergent dissi¬ 
pation results as long as they include a large enough 
number of low Ethr daughter and granddaughter modes. 
Moreover, it is not necessary to include collective sets of 
daughters and granddaughters; although they can mod¬ 
ify the dynamics somewhat, the low Ethr three-mode sets 
ultimately model the total dissipation well. 

Further numerical experiments reveal that a network 
consisting of one parent, its ~ 20 lowest Ethr daughters, 
and ~1500 low Ethr granddaughters is sufficiently large 
that it yields convergent results while still allowing us to 
efficient explore a range of Mp and P (see caption of 
Figure 1^. We use this as our “reference network” when 
computing {te) as a function of Mp and P in §[^ That 
such a network is sufficiently large can be gleaned from 
the left panel of Figure]^ which shows that the reference 
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Figure 7. Mode energy Ea (left panels) and effective number of modes A^eff (right panels) as a function of time for reference networks 
with Mp = Mj. (top to bottom) P ~ 1 day on resonance, off resonance, P ~ 3 day on resonance, off resonance. For clarity, on the left we 
plot only one out of every 25 granddaughter modes included in the simulation. 


network is very similar to A^gen = 3 or 4 networks 
with iVmodes >10^. 

We demonstrate this further in Figure which shows 
the mode energy and the effective number of modes par¬ 
ticipating in the dissipation as a function of time 
for four different orbital periods. We estimate Ngg by 
computing 


where 


iVeff = exp \ - ^Pa Inp, 


Pa = 


laQaC 
Ea laQala ' 


(23) 


(24) 


This statistic is relate d to the Shannon e ntropy and is 
similar to one used in Brink et al. (2005). If all modes 
contribute equally to the dissipation, then = l/V„iodcs 
for each mode and Neg = N^oges- We see in Figure 0 
that the dynamics are complicated, with many excited 


modes. We also see that the mode energy and N^g in¬ 
crease at shorter P and near linear resonances. How¬ 
ever, Neg Vmodes while both the peak energy of the 
modes and N^g remain nearly constant after ^100,000 
orbits, indicating that our reference networks are suffi¬ 
ciently large. The behavior of the networks as a whole 
is what is important here, rather than the dynamics of 
any individual mode. Individual modes fluctuate signifi¬ 
cantly but the overall dynamics do not change over long 
timescales. 

We now describe our procedure for calculating {te) as 
a function of Mp and P, the results of which we present 
in §[^ For each reference network run at a given (Mp, P) 
point, we simulate at least 5 x 10® orbits in order to al¬ 
low transient effects from initial conditions to die awayj^ 

^ This is true of all networks in Figurel^as well, with the excep¬ 
tion of the largest Ngens = 4 network ana the collective network, 
where we were computationally limited to shorter (but still con¬ 
vergent) integrations. 
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Figure 8. Dependence of {te) on P for fixed values of Mp (left panel) and the dependence of {te) on Mp for fixed values of P (right 
panel). The solid lines show our analytic fitting formula, which begins to break down for Mp < 0.5Mj. 


Using Equation (11), we then compute the average E* 
over the last ~ 5 x 10^ orbits of the integration. We 
do this in order to average over the rapid fluctuations 
in dissipation that characterize the nonlinear equilibria, 
which correspond to r.m.s. variations at roughly a 10% 
level. We also find that E* depends slightly (factor of 
« 2) on how close the parent is to a linear resonance 
peak (see the iVgen = 3 results in Figure]^. In order to 
compute the average dissipation rate, we must therefore 
average the results over several resonance peaks. We do 
this by performing 21 separate integration runs, each at 
a slightly different orbital period {SP ^ P) chosen such 
that the runs span three resonance peaks. We then cal¬ 
culate the average E* weighted by the amount of time 
spent at that orbital period (Appendix]^ and use this 
{Et) to estimate {te) via Equation (13). 


5. RESULTS 
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In this section we present the results of integrating the 
coupled amplitude Equation using the procedure and 
reference networks described m § |4.6[ Figures and 
show (te) as a function of Mp and P. We find that 
{te) depends strongly on P and only mildly on Mp. Our 
numerical results are well approximated by the fit 


over the range 0.5 < Mp/Mj < 3 and < 4 days. 
This matches our numerical results to within a factor 
of ~ 2 over this range, which is comparable to system¬ 
atic modeling uncertainties due to differences in the net¬ 
work structure of larg e A^gens = 3,4 networks (see Figure 
[^. By Equation ([14), this corresponds to a stellar tidal 
quality factor 


Q(. = 3.0 X 10^ 


(Mp\ 

\M,) 



(26) 


Thus, for Mp > 0.5Mj and P < 2 day we find that (te) is 
small compared to the main-sequence lifetime of a solar- 
type star. For Mp < 0.5Mj, we find that although (te) 
increases significantly, it can still be small at small P. 


Figure 9. Decay time {te) on the {Mp, P) plane. The size 
of the marker is proportional to logj^Q (Qt), with corresponding 
color bar for {te) in years. Each sample poi nt is labeled with 
l ogin f ( tk i/Gvr). The shaded region represents |Barker & Ogilvie| 
| |2011[ l’s prediction for when linearly resonant modes break. 

For example, for Mp — 0.1 Mj and P — 1 day we find 
(te) — 800 Myr. At sufficiently small Mp and/or large 
P, the non-linear effects “turn off” and (te) collapses 
onto the linear result. This is particularly evident for 
Mp = O.lMj in Figure]^ We also illustrate this effect 
in Figure [l0[ which shows how (te) depends on P when 
Mp = Mj for different numbers of generations. We see 
that by P ~ 10 days, {te) is very long and close to the 
linear prediction (the Ng^n = 1 line). 

5.1. Implications for a few known systems 

Based on the Exoplanet Orbit Database 
(http://www.exoplanets.org), there are currently 
11 known planets orbiting approximately solar-type 
stars (M = 1.0 ± O.IMq and Tgff ~ 5500 K) with 
decay times {te) < 1 Gyr according to our results. Of 
these, 7 have expected decay times {te) ^ 0.3 Gyr; in 
order of increasing P, they are WASP-19b, TrES-3b, 
HAT-P-36b, WASP-77Ab, WASP-4b, WASP-36b, and 
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Table 1 

Orbital parameters for known systems and a summary of predictions for the orbital decay timescales and change in orbital parameters 

after 10 years 



WASP-19b 

HAT-P-36b 

WASP-36b 

CoRoT-2b 

M[Mq] 

0.930 ±0.02 

1.022 ± 0.049 

1.020 ±0.032 

0.970 ±0.06 

rIRq] 

0.990 ±0.02 

1.096 ±0.056 

0.943 ±0.019 

0.902 ±0.018 

Mp sin i [Mj] 

1.114 ± 0.039 

1.83 ±0.1 

2.255 ± 0.089 

3.27 ± 0.171 

P[sec] 

68155.776 ± 0.026 

114682.78 ±0.26 

132828.36 ± 0.23 

150594.64 ± 0.86 

eccentricity 

0.0046t°;°°« 

0.063 ± 0.032 

0 

U.U140_Q QQ7g 

age[Gyr] 

io.2l®;° 

null 

null 

Q y+3.2 

-2.7 

{te) = E,rb/{dtEA [Myr] 

9.2 ± 0.128 

205 ±3.9 

454 ± 15.7 

623 ± 27.6 

min TE [Myr] 

6.3 

84.7 

214 

241 

max Te [Myr] 

12.4 

311. 

853 

1150 

AP = (3P/2r£;)At [ms] 

110. ± 1.5 

8.39 ± 0.16 

4.4 ±0.15 

3.6 ±0.16 

min AP[ms] 

82 

5.5 

2.3 

1.96 

max AP[ms] 

161. 

20. 

9.3 

9.38 

Tahift = (3/4TB)(At)2 [sec] 

257 ± 3.6 

11.5 ±0.22 

5.2 ±0.18 

3.80 ±0.17 

min Tahift [sec] 

191 

7.6 

2.8 

2.1 

max Tahift [sec] 

375 

27.9 

11. 

9.8 



P[day] 


Figure 10. Decay time {te) as a function of P for different 
numbers of generations and Mp = Mj. (blue squares) 10 most 
resonant parents, (red circles) parents+daughters. (green trian¬ 
gles) parents+daughters-|-granddaughters. (dashed lines) analytic 
estimates for Nggns = 1,2 (Appendix]^. 


WASP-46b. Figure shows these planets on the 
/ith 


MpSini-P plane, witli (r^;) labeled for each system 
and a contour of (te) superimposed. These planets all 
have Mp sin i > Mj, P < 2.0 days, and eccentricities 
consistent with or very close to zero. Since these are 
all transiting systems, Mp sin i ~ Mp and the reported 
errors in the measured mass are typically < O.lMj. 

We note that two of the planets (CoRoT-2b and 
CoRoT-18b) have masses MpSini > 3Mj. This sug¬ 
gests that they are in the strongly nonlinear reg ime where 


the parent wave breaks within the stellar core (Barker & 
Ogilvie 2010 2011[ [Barker 2011). We discuss how our 


results compare to the strongly nonlinear simulations of 
Barker & Ogilvie in § [^ 

Of the 11 planets with (te) < 1 Gyr, there are five 
for which studies report at least some constraint on the 
age of the system. In three of these, the age uncertainties 
are sufficiently large that the systems might be relatively 
young, i.e., ~ 1 Gyr (WASP-64b, WASP-5b, CoRoT-2b: 
1 .2tJ;7, 2.712^7, respectively). However, WASP- 

4b and WASP-19b appear to be older systems: 7.0 ± 


3.5 


3.0 


g2.5 

I 2.0 


1.5 


1.0 


Figure 11. Period-mass distribution (Mp sin 2 versus orbital pe¬ 
riod P) of known extrasolar planets orbiting solar-type stars with 
{'^e) $ 1 Gyr based on our res ults . Each planet is labeled by our 
(te) fitting formula (Equation [25|l , and the color represents con¬ 
tours of {te)- (Data taken from 1ittp://www.exoplanets.org). 
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2.9 and 10.2)')|'gGyr, respectively. Assuming that the 
planets arrived close to their current orbits when their 
host stars first formed, such old stellar ages seem to be 
in tension with the small (te) we predict, especially in 
the case of WASP-19b. If our results are correct, then 
perhaps these planets were scattered into their current 
orbits well after the stars formed or they just happened 
to initially reside at separations with decay timescales 
very close to their current ages. 


Several recent papers consider the prospects for the di¬ 
rect detection of orbital decay of individual planets by 
measuring transit timing variation s (TTVs) over long 


time 

baselines (At > 5 year, see Gandolfi, D. et al.| 

20151 

IBirkbv et al.||2014| jValsecchi & Rasio||2U14MWat-| 

son & 

Marsn|2010). In order to evaluate this possibility. 


we simulated four known systems spanning a variety of 
companion masses and orbital periods (but each with a 
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Figure 12. Cumulative probability distributions for observable changes as a function of elapsed time, (blue) WASP-19b. (red) HAT- 
P-36b. (green) WASP-36b. (purple) CoRoT-2b. Solid lines represent approximate fits to the simulation results, which are represented by 
filled circles. 


solar-type hos10, calculating their tidally induced TTV 
(Tshift) and change in orbital period (AP) as a function 
of At . We compute these according to (see Birkby et al. 
2014 for a derivation) 


AP « PAt = = — At. 

2 Porb 2r£; 


(27a) 

(27b) 


I9b (jHebb et al. 

2010 Hellier et al. 2011 

Mortier 

et al. 1^01(11), HA'i'- 

L'-3(jb (Bakos et ai. 2U12D, WAyP-^tib 

([Smith et al.| 

2012 

), and (JoKo'i'-2b (Alonso et al. 2008 

Gillon et al.||2010D. 

In order to calculate the orbital decay rate of these 


systems, we simulate a small range of orbital periods 
centered on each system’s measured period. We then 
compute the time-averaged decay rates following the 
procedure described in § |4.6| We do this in order to 
mitigate any differences between the resonances of our 
stellar models and the actual resonances of the stellar 
hosts. Furthermore, this allows us to compute a mini¬ 
mum and maximum expected decay rate, corresponding 
to the troughs and peaks of the resonances, respectively. 

Table lists (te) as well as the minimum and max¬ 
imum te- The {te) of the four systems ranges from 
about 10 Myr (WASP-19b) to 600 Myr (CoRoT-2b), 
while the minimum (maximum) te is approximately two 
times smaller (larger). WASP-19b has by far the short¬ 
est decay time owing to its extremely short orbital period 
(18.9 hours). 

Table also lists the systems’ average, minimum, and 
maximum AP and Tghift after ten years of evolution. 
These provide an estimate of the magnitude of the tidally 
induced deviations we would expect to observe from these 
systems over the next ten years. 

We quantify these effects further in Figure which 
shows the cumulative probability of observing tide- 

® This require ment is why we do n ot consider WASP-18b, which 
was analyzed in [Birkby et al.| ||2014|. 


induced deviations as a function of time. We choose a 
detection threshold of (AP)^jjj, = 0.1 sec and (rshift)thr = 
60 sec based on the measurement errors of P and the ex- 


pected uncertai nties in TTVs (Gillon et al 


1200 ^ 

lie At 


Watson 

through 


& Marsh 2010); different choices will sea. 

Equation (|27|r 

We find that TshUt should always produce a detection 
faster than AP. This is because Tghift is a cumulative 
effect that builds up throughout the orbital decay. Ac¬ 
cording to our results, WASP-19b should produce a de¬ 
tectable Tghift in the very near future, with a « 50% 
chance of observin g a deviation now given the current 
four year baseline (Hellier et al. 2011) and a high likeli¬ 
hood of detection after only two more years. It will take 
considerably longer before detections are possible in the 
other three systems. 

We note that even if, for some reason, our calculations 
overestimate the dissipation ra te b y an order of magni¬ 
tude, the Tshift curves in Figure [T^ would only b e shifted 
to the rig h t by a factor of vTO ~ 3. Finally, as 


Watson 


& Marsh (2010) point out, the Applegate effect could 
produce AP and Tshift values that are comparable to the 
tidally induced values and distinguishing the two may 
not be simple. 


5.2. Comparison with previous estimates of nonlinear 
tidal dissipation 

Previous studies that attempt to estimate the nonlin¬ 
ear dis sipation rate of dynamical tides in close binaries 
includ e Kumar & Goodman (|1996 ) and Barker & Ogilvie 
(2011). I'hey both argue that an upper bound to the dis- 


sipation rate is approximately given by the product of the 
parent’s linear energy and the three-mode growth rate of 
the fastest growing daughter pair: 


T* ^ PsmdTlin. 


(28) 


This estimate does not account for the continuous lin¬ 
ear driving of the parent. Instead, the parent wave is 
initialized with an energy equal to Pun but is otherwise 
undriven, and the problem reduces to determining the 
amount of time it takes for daughters to dissipate that 
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initial energy. Althou gh this is appropriate for the tidal 
capture problem that Kumar & Goodman (1996) con¬ 
sider (because the bihary is on a highly jccentric or¬ 
bit and the parent is only driven strongly during the 
brief pericenter passage), in our analysis the orbit is cir¬ 
cular and the parent is a contin uous ly driven standing 
wave. The estimate of Equation (28) also assumes that 
the mode dynamics are dominated by the single, fastest 
growing daughter pair even though there may be many 
modes participating in the interactions. 

By Equation (17), we find that the fastest growing 


daughters have a growth rate 


Esmd — 0.6 


- 11/6 


(-) yr- 

mJ Vday; ^ 


(29) 


and by Equation (15) 




erg s 


( 30 ) 

For comparison, the fit to o ur numerical simulations 
yields, by Equations (25) and (12), 


E* d. 3.5 X 10^9 

Mj) Vday 


-7.4 


erg s 


-1 


(31) 


Thus, while the two have nearly identical P scalings, 
the E* from our simulations is ~ 15 times larger than 
TamdEiin for Mp ~ Mj. This factor of 15 difference 
can be seen in the Q'^ estimates. In particular, we find 
Q* « 3 X 10® at P = 1 day for Mp > 0.5Mj. By con¬ 
trast, Barker & Ogilvie (2010) argue that Q* > 5 x 10® 


for systems below the wave breaking limit {Mp < 3Mj) 
based on their assumption that P* < EamdPiin in the 
weakly nonlinear regime. 

We suspect that the discrepancy is largely due to the 
assumption in Equation (28) that only the single fastest 
growing daughter pair is important. In Figure we 
demonstrate that this is not the case. We show tEe in¬ 
dividual and cumulative contribution to P» of modes in 
our reference network (which consists of ~ 1500 modes). 
We find that there are several daughter modes that con¬ 
tribute substantial amounts of dissipation, not just a sin¬ 
gle dominant daughter pair. Figure also shows that, 
in sum, the granddaughters are the dominant source of 
dissipation in the network. 


6. SUMMARY AND DISCUSSION 

We present a first principles calculation of the satura¬ 
tion of nonlinear interactions between g-modes excited 
within the cores of solar-type hosts by planetary com¬ 
panions. Using a WKB approximation for high-order, 
adiabatic g-modes and analytic approximations to their 
coupling coefficients detailed in WAQB, we systemati¬ 
cally investigate the number of modes and types of cou¬ 
plings that are dynamically relevant. We determine the 
minimum mode network size and structure that yields to¬ 
tal dissipation rates consistent with those of much larger 
networks (to within a factor of « 2). This minimum net¬ 
work is sufficiently nimble that we can efficiently explore 
broad swaths of the (Mp, P)-plane. We find that weakly 
nonlinear interactions are energetically important over 



Figure 13. Energy dissipation rate Ea for each mode in our ref¬ 
erence network, which consists of one parent (blue square), ~ 20 
daughters (red circles) and ~ 1500 granddaughters (green trian¬ 
gles). Modes are ordered by Ea and here we take Mp = Mj and 
P ~ 3 days. The blue line is the cumulative distribution of Ea‘, we 
see that granddaughter modes are responsible for the majority of 
the dissipation. 


large portions of this plane, including regions occupied 
by known exoplanetary systems. In these regions, the 
orbital decay time (te) and stellar tidal quality factor 
Q' follow simple power law relations (Equations and 


26). 


We find that the orbital decay of a number of observed 
hot Jupiters should occur on timescales much shorter 
than the main sequence lifetime of their host star. Such 
rapid orbital decay could explai n the observed paucity 
of giant planets with P < 2 day s (McQuillan et al.]|2013 


see also |Winn fc Fabrycky|20T5] for a recent review of the 
observations). The short decay times would also induce 
TTVs that may be observable with current technology 
(especially that of WASP-19b). Precision photometry of 
individual systems may thus provide a new handle on 
tidal interactions within the next few years. 

Our calculation comes with some caveats. First, al¬ 
though our reference network yields dissipation results 
that are very similar to those of the largest networks we 
investigate (which have > 10 times more modes than the 
reference network), there is still a possibility that the 
dynamics will change upon the addition of even more 
modes. Second, our calculation assumes that the modes 
are all global standing waves. However, this prescription 
may break down if the amplitudes of the modes change 
on timescales shorter than the group travel times be¬ 
tween their inner and outer turning points. Moreover, 
although the parent mode is below the wave breaking 
threshold (when not too close to a linear resonance), the 
daughter and granddaughter modes may not be. In Ap¬ 
pendix]^ we show that the threshold amplitude of the 
three-mode parametric instability is much smaller than 
the wave breaking threshold and that both have the same 
frequency scaling. This may mean that further gener¬ 
ations will be excited before the daughter and grand¬ 
daughter modes break. Nonetheless, this issue deserves 
further investigation. Finally, we do not account for pos¬ 
sible changes to the stellar structure due to the transfer 
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of energy and angular momentum from the sea of ex¬ 
cited waves. Further work is needed in order to deter¬ 
mine the extent to which stellar spin-up, heating, and/or 
evolution affect background properties such as the star’s 
stra tification and thereby t he wave interaction dynamics 
(see Barker fc Ogilvie|2010[ for a discussion of this issue). 

Our study focuses on wave interactions in the weakly 
nonlinear regime. For solar type stars, this corresponds 
to planetary masses Mp < 3.6Mj(P/l day)“°-^; above 
this mass, the parent wave breaks as it approaches the 
stellar center and the system is therefore in the strongly 


nonli near regime (Barker & Ogilvie 2010 2011 Barker 
20111. In the weakly nonlinear regime the parent is 
a global standing wave while in the strongly nonlinear 
regime the parent is more appropriately treated as a trav¬ 
eling wave; it does not reflect upon reaching the stellar 
center. Barker & Ogilvie study the fate of such a strongly 
nonlinear traveling wave with numerical simulations us¬ 
ing a Boussinesq-type model. Because our calculation 
studies a different hydrodynamic regime, a direct com¬ 
parison with their results is not possible. Nonetheless, 
one might expect the two to roughly agree near the re¬ 
gion that marks the transition from weakly nonli near to 


strongly nonlinear (i.e., near Mp ~ 3Mj). Indeed,[Barker 
(20111 finds Q* « 9 x 10^(P/day)^ ® for waves that break 
while we find Q), « 5 x 10®(P/day)^ '^ for Mp ~ 3Mj. We 
explore some of the similarities between the two regimes 
further in Appendix [F} 

We find Q/ « 3 x 10® at P = 1 day for Mp > 0.5Mj. 


This appears to conflict with the estimate in [Barker fc 
Ogilvie| (|2010), who argue that Q* increases rapidly to 


10 ® for systems below the wave breaking limit 
(Mp < 3Mj). They do not attempt to calculate the 
saturation of the nonlinear parametric instabilities as in 
our study but instead base their estimate on stab ility 
analysis scaling arguments. As we explain in § |5.2[ the 
issue might be that their estimate neglects the continu¬ 


ous driving of the parent and does not account for the 
complicated multi-mode dynamics that we find are im¬ 
portant. Interestingly, we do see a steep increase in {te), 
although at much lower Mp. 

In order to be consistent with the observed distribution 


of exoplanets, Penev et al . 


]2012) find that > 10^. 


However, as [Birkby et al. (29141) note, their analysis is 


for one specihc set of initial conditions with some ideal¬ 
ized assumptions about the chances of a planet candidate 
being confirmed by follow-up. They also assume gas disk 
migration and, as [Penev et al.j point out, their result may 
not be valid for other giant planet migration mechanisms 
such as dynamical scattering. If gas migration is the 
dominant mechanism that creates hot Jupiters, then our 
results suggest that hnding these systems at P < 2 days 
should be extremely rare. However, if scattering popu¬ 
lates short period orbits at random times after a system’s 
formation, then a low may not necessarily conflict 
with the observed population of hot Jupiters orbiting 
^ Gyr old hosts. 

Our study only considers solar-type hosts even though 
the observed population of hot Jupiters includes a wide 
variety of host types. Since the linear and nonlinear ex¬ 
citation of waves by the tide is sensitive to the detailed 
structure of the star, it is not clear how our results might 
depend on stellar type. Extending the analysis to non¬ 
solar type hosts would therefore allow us to more fully 
assess the prospects for measuring tide-induced orbital 
decay of individual hot Jupiter systems. 
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APPENDIX 


A. COMPUTATION OF TIME-AVERAGES 


Quantities such as the energy dissipation rate E* depend on how close the system happens to be to the densely 
spaced linear resonance peaks (the frequency spacing is |Aa|/H « 10“^(P/ day)“^). Because we are mostly interested 
in time-averaged statistics, at each (Mp, P) point, we carry out 21 distinct simulat ions , each separated slightly in 
orbital period with a spacing chosen such that they span three resonance peaks (see §14.Gland Figure]^. We compute 
the time-averaged statistic of a quantity X by weighting each sample by the amount of time spent at that period 


fdtx fdpp-^x 
Jdt f dPP-i 


(Al) 


where P = dP/dt is the rate at which the period changes due to tidal dissipation. We compute P using an energy- 
balance argument. We expect the time rate-of-change of the orbital energy (Porb), the tidal interaction-energy (Bint), 
the rotational energy of the synchronized companion (Prot), and the energy stored in the modes (Pmodes) to balance 
with the energy lost through dissipation 


d 

dt 


(Aorb d~ Bint d~ Erot d~ Emodes) 


i 


(A2) 
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from which we can extract the time rate-of-change of the orbital period via 



'' dPorb 

, dPint , d 

. dP 

+ dp +“ 

-H 

'dPorb 
. dP 

) 

' i 

p 




-1 


dP 


dP / 




Porb 


(A3) 

(A4) 

(AS) 


where we have noted that |dPorb/dP| ^ |dPint/dP|, |dProt/dP|, |dP„iodes/dP|. This is because Porb is much larger 
than any of the other energy scales, so even small relative changes in Porb dominate over the other terms. This gives 


(A) = 


/dP 


P 

Earh 


-1 / \ -1 
P 

So.b 


/dP 




Ea- 


J2liEt 


E 

P 




A 


(A6) 


We use this procedure to calculate the time-averaged P» in the neighborhoods of each orbital period and {te). 

B. THREE-MODE NON-LINEAR EQUILIBRIUM 

Here we briefly review the non-linear equilibrium for three-mode systems. The calculation is similar to that of 
Appendix D of WAQB except that here we provide more detail about the phase relations amongst the modes. We 
begin with the equations of motion (Equation and introduce the change of coordinates x = yielding 


dtxp + iiAp+^p)xp = 


(Bla) 

(Bib) 

(Blc) 


We can cancel all time dependence in these equations by demanding 

rriafl = UJa — Aa = Ap + A..^ — LOp — LO^ (B2) 

and assuming that i9t —>■ 0 in order to explicitly seek time independent solutions. Manipulating the two daughter 
equations yields 


[lAp + 'Yp)xp{-iA^ + 7^) = 2iuipkap^x*^ + Jj)xj)* = Aujpuj^kip^XaX*^xp 


which implies 


Ap-f^ = A^jp, 


ApA^ + 7/377 — 


where we write x = Ae*°. We then have 


^2 ^ A^A^ + 7^37^ 


7/377 




1 -L 


Ap + A^ 

7/3+77 


and we recover the parent instability threshold energy Pthr = (Equation [^. The daughter equations yield 


{iAp + 'Yp)xp _ {iA^ + 'y^)x^ 


UJpXj 


OJjXp 


Af- 


77^/3 

7/3W^’ 


which gives 
or equivalently 


{iAp+-fp)Ap = 2itupkc,p^Ac.A^e-^^^-+^^+^X^ 


Ar- 


7/3 — = 2a;/3fca/37AQsin(5, 


Ar 

Ap— — 2ujpkctP'yA(^ cos(5, 


(B3) 

(B4) 

(B5) 

(B6) 

(B7) 

(B8) 


where S = Sa + Sp + S^. We can now use the parent equation to determine the parent phase Sa and the product of the 
daughter amplitudes 

{iAa + 'ya)Aa = iu)aUae~^^°‘ + 2ioJakap’yApA^e~^^. (B9) 
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After some manipulation, we find 






oc^a.^'y 


{Aa COS ^ + 7 q, sin i5) ± W (cos 5 + 7a sin i5)^ + 


A 2 


A 2 


a;2[72 

rv rv 


Al 


-Al 


■7a 


(BIO) 


The choice of sign depends on the sign of uJanap-y and is determined by the requirement that the daughter amplitudes 
be positive. Finally, by Equation (15), we see that the instability condition is 


^'oJAg _ Alin »2 _ Athr 

A2 +7^ “ Ao ^ Ao ■ 


(Bll) 


Note that we can solve for the parent’s phase 5a and the sum of all the mode phases (5, but we cannot break the 
degeneracy between the daughters’ phases. This is observed numerically, and carries information about the initial 
conditions. 


C. TWO DAUGHTERS, N PARENTS 

If we linearize around the linear-equilibrium solution, the equations of motion for the daughter modes become 

dtqp + (iuj/s + 7/3) = 2iujpq* ^ Kp/B^qp (Cl) 

p G parents 

for daughter /3 and the equivalent equation with the exchange j P for daughter 7. We can analyze this system as 
if there is a single parent with complex amplitude 

Kq = 'y ( Kpfj^qp. (C2) 

p G parents 


We note the possibility for parent modes to interfere with one another when driving daughter modes, possibly 
rendering daughters stable under multi-parent driving when they were unstable to any individual parent. Most 
notably, if the parents are nearly regularly spaced in frequency and driven at the midpoint between their resonance 
peaks, there can be strong destructive interference. This is because each parent is paired with a partner on the opposite 
side of the driving frequency, and each pair consists of parents oscillating with nearly opposite phases. This narrow 
“trap” in the resonance troughs is readily apparent at orbital periods above 4 days for a solar-type host of a Jupiter 
mass companion. However, we did not observe significant “trapping” below ~ 4 day orbital periods, where we focus 
our attention for this study. This may be due to the asymmetric spacing of resonances, which will destroy this near 
perfect cancellation, or due to the amplitudes being large enough to overcome any cancellation that was present. In 
the hot Jupiter context, this issue is probably only of theoretical interest since the orbital evolution time scales are 
> 10^^ yr for P > 3 day, even for a 3Mj companion. 


D. DETAILS OF COLLECTIVE SET SELECTION ALGORITHM 


One can easily think of more complicated collective sets than what is described in § |3.2[ We analyze several of these 
systems in Appendix]^ In order to detect and include the diverse set of collective systems, we implement a broad 
search through parameter space. We begin with a “seed” triple in parameter space, typically taken to be a minima 
of Athr- We then expand the set of included modes in (n, Z)-space around these seeds, choosing new modes from the 
border of the included set. For these border modes, we compute the three-mode Athr for all possible couplings between 
that border mode and the interior modes. We then sort these Athn and divide each by the number of couplings that 
produce Athr less than or equal to the current value. We take the minimum ratio and call it the “collective Athr-” This 
approximates the s caling with N predicted in Appendix and incorporates the decoupling of large detuning modes 
discussed in § |E.4| Border modes are added in order of increasing collective Athr, and these thresholds are updated 
each time a mode is added to the network. If the detuning increases the three-mode Athr faster than the number of 
modes included, then small sets with low detuning will naturally be chosen. However, if the detuning increases Athr 
more slowly than the number of modes included, then the collective Athr will decrease with the addition of more modes 
and the algorithm will select a set of collectively unstable daughters. 

We typically find that a minimum number of daughters is needed before the scaling with N dominates over the 
increase in Athr- Depending on Ahn, these collective sets can grow to sev eral thousand modes. Although each mode 
can only directly couple to a relatively small number of other modes (see § |3.2| and WAQB), we find that many sma ller 
sets overlap and are thereby strung together to create larger networks. We discuss some of this behavior in_§_ E.2 

This algorithm scales poorly with the number of modes included (0[A^]). Eurthermore, as we describe in we find 
that we can accurately model the total dissipation within the star using only three-mode pairs, rather than collective 
sets. This, coupled with the fact that large collective networks are expensive to integrate, is the re ason we ch oose 
three-mode networks with many couplings and generations as our reference networks discussed in §§ |4.6[ [5| and |5.1 
We note, however, that collective sets may be important if one is interested in accurately modeling the dynamics o: 
any particular mode, rather than the network as a whole. 
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E. COLLECTIVE SET STABILITY 

In this appendix we analyze collective instabilities, i.e., sets of daughter modes that display rapid growth rates due to 
their mutual inter-coupling. Although in our simulations we find that they do not contribute significantly to the total 
tidal dissipation in hot Jupiter systems, for completeness we present here derivations of different stability thresholds 
for different types of collective instabilities. Our mode selection algorithm (Appendix]^ finds complicated collective 
sets that contain these types of coupling topologies. 


E.l. Single collective set 

We first consider the stability of a single collective set. Since we are interested in the stability of linear solutions, we 
can assume that the parent is at a fixed amplitude 


go = (El) 

The equation of motion of each daughter is then 

Qa + (iuJa + la)qa = (E2) 

/3 

= -I- 2i0Ja ^ K,oal3Aoe'^'‘^^*~^°\p. (E3) 


Defining a new set of variables q = xe we can re-write the daughter equations as 

A T* —2Aa)i— 

-\-0J^ — ^/3 


Xot “1“ “1“ j3^ 


-\-‘2iUJcy E 




2(5o 






where in the last step we demanded that the time dependence cancels 

^ ‘^/3 — — ^/3 = 0 V {o!, /3}. 

Analyzing this as an eigenvalue problem, separate x into real and imaginary parts x = R 


(E4) 

(E5) 

(E6) 


(E7) 


Ra 


H” ^oi-^o^ocyoL sin So Aq, “h cCq Aq/^qqq cos Sq 


Ra 

.4. 


Aq “h CCq AqAvqqq cos Sq Oi AqA^^qq Sln Sq 


.Ic^. 


+ E 

/3/a 


2uJaAol^oafS sin So 

2uJaAo^oal3 COS So 


2iuq. AqNock/? COS Sq 
2 iOoiAoKoai 3 sin So 


If we assume [Raila] oc V a, then this equation becomes 


Rp 

Ip. ■ 

(E8) 


(^Q H“ s) ~h Aq^^oqq sin Sq 

Aq ~h ^CkAq/^oqq cos Sq 

Ra 

Aq “h WqAo^Yoqq cos Sq 

-{lot + 5) - UJotAoKootOL sin So_ 



2ujaAoKoap sin So 
2 u}oiAoKoal 3 cos So 


2wq Aq/^oq^ cos Sq 


Rp 

- 2 u;o,Aonoa 0 sin So^ 


[lp\ 


(E9) 


This is an eigenvalue problem for a large matrix and the general decomposition is difficult. However, the matrix 
can be made almost symmetric and if we make several approximations the problem becomes analytically tractable. 
Specifically, if we assume 


Wq , = w 

7a = 7 


V a. 


^oaP — OL ^ {5^ 


(ElO) 


then we can define 


Ms 


~ (7 + s) + cjAoKs sin So 
—A -I- ujAoKs cos So 


A -I- LvAoKs cos So 
— (7 -I- s) — ijjAoHs sindoj ’ 


Mi = 


2a;AoKsin(5o 
2 tLiAo«; cos So 


2ujAoK, cos So 

— 2 L>jAoKS\nSo ’ 


(Eli) 


where Aq, = A V a since Wq = a; V a. Writing this as a single matrix and requiring non-trivial mode amplitudes, we 
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obtain 


0 = det 


Ms 

Mi 
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Mi 

Ms 
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Mi 

Ms 
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■ Mi 

Mi 

Mi 

Mi 

Mi 

Mi ■ 

■ Ms 

Mi 

Mi 

Mi 

Mi 

Mi ■ 

■ Mi 

Ms 


= det \Ms - det \Ms + {N - 1)M/|, 


(E12) 


We have N — 1 repeated pairs of roots and one additional pair. The eigenvalues can be easily computed from 

det \Ms — Mj\ = (7 + s)^ + — 2k)‘^ = 0 

s = —7 ± — Ks)2 — 

and 

det \Ms + (A - l)Mi\ = (7 + s)2 + A2 - uj^Al(2{N - 1)k + Ksf = 0 
=> s = —7 ± A^{ 2 {N — 1 )k + KsY — A^. 

In particular, we are interested in the values of Ao for which ]R{s} —>■ 0. These are 


- 

^thr “ 


■A^ 


77 


and 


A^ - 

"^thr — 


4uj^(k—^Ks)^ 4ujuj(k — ^Ks)^ 
72 + A 2 


(A + A)2 

(7 + 7)2 


4uj^((N — 1 )k + \N — 1J 4ujuj(k + 


77 


2(A-1) 


1 + 


(A + A) 
(7 + 7 )^ 
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(E13) 


(E14) 


(E15) 


(E16) 


respectively. We see that there are A — 1 modes that resemble “standard” three-mode instabilities and one collective 
eigenvalue, with an amplitude threshold su ppress ed by a factor of A — 1. 

Because of the assumptions in Equation (ElO), the actual value of ^thr will differ somewhat from this expression. 


Nonetheless, we expect it to generalize to the requirement that 


(N-lfA^ > 


7172 


4wia;2Ko 


(Ai 


(71+72)^ _ 


V modes 1, 2 S collective set of A modes. 


(E17) 


where Ai -|- A 2 = + oji + 0 ^ 2 - 


E. 2 . Overlapping collective modes stability 

We now consider a coupling topology where there are three types of modes. The A modes are coupled to other A 
modes and to C modes. B modes are coupled to other B modes and to C modes. C modes are coupled to all other 
modes. Furthermore, we assume that all A, B, and C modes are coupled to the same parent modes, which we treat 
as a single parent even though multiple parents may be acting (see Appendix [^. 

The associated eigenvalue problem yields the following characteristic equation 
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(E18) 
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We again note the high degree of symmetry, which allows us reduce the determinant to 

0 = (det \Mi - (det |m| - Mf (det \M'§ - 

X det \Mi + {Na - 1)M/| det |M| + {Nb - l)Mf | 

xdetjMf+ (fVc-l)Mf(E19) 

where 

Z = NaMf {M§ + {Na - Mf + NbM^ {M§ + {Nb - Mf. (E20) 

We recognize this as Na — I independent A eigenvalues, Nb — 1 independent B eigenvalues, Nc — 1 independent C 

eigenvalues, one eigenvalue corresponding to the collective modes without the coupling to C modes for each of the A 

and B modes, and a collective set for the C modes with a modification due to the couplings to the A and B modes 
(through Z). We further note that when Nc —>■ 0, the eigenvalues reduce to two separate collective sets, as expected. 

The interesting eigenvalue is due to the interaction between the C modes’ collective set and the couplings to A 
and B modes. If we assume that all mode parameters are the same for all sets of modes, and further assume that 
Na = Nb = Nc-, we can make analytic progress on this determinant, and obtain 

(7 + s)" + - uj^Al {{ks + 2{N - l)kf + 8N^k^) = 

and the threshold amplitude 

a2 _ 7^ + 

“ 4a;2 {3k^N^ + k{ks - 2k)N + k{k - k,,) + k^j4 
7" + A 2 \ 1 /72 + A 2 

4a;2fc2 ) Nf+N^+N^ 

where we assumed the limit of large N. We note that this is very similar to the case of a single collective set, except 
A2 —)• N^ + N^ + Nc ■ If we stitch together many separate collective sets by overlapping them, we only expect the 
effective number of modes to sum in quadrature. This was tested numerically by taking the determinant without 
assuming equal numbers of modes, and found to be in reasonable agreement with this scaling. 



= 0 


(E21) 


) 


(E22a) 

(E22b) 


E.3. Non- “self coupled” collective sets 

Appendix |E.1| and E.2 considered self-coupled modes. However, the vast majority of couplings will be between 
modes that do not support self-coupled daughters. For example, if the parent azimuthal order m is odd, then the 
daughter modes must have different m numbers. If we consider two sets of modes, one with N daughters and one 
with n daughters, we can define 2x2 sub-matrices similar to Appendix |E.1| for each group of modes. This means we 


will also find collective sets with characteristic equations like the following, with capital letters corresponding to the 
A-mode set and lower case letters corresponding to the n-mode set 


Ms 

0 •• 

■ 0 

0 

Mi 

Ml ■■ 

• Mi 

Mi 

0 

Ms ■■ 

• 0 

0 

Mi 

Mi ■■ 

■ Mi 

Mi 

0 

0 •• 

• Ms 

0 

Mi 

•• 

■ Mi 

Mi 

0 

0 •• 
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0 • • 
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M, ■■ 

■ M, 
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0 • • 

■ 0 

Ms 


where this is an {N -|- n) x {N -\- n) matrix. We can simplify this to only a 4 x 4 determinant 


(E23) 


0 = (det |Ms|)'^ ^(det|Ms|)” ^ det 


Ms 

NM, 


nMj 
Ms ’ 


(E24) 


which looks like a set of independent eigenmodes and a 4 x 4 determinant for the collective modes. In general, that 4x4 
determinant must be solved numerically. However, if we again assume identical mode parameters and that N = n, we 
see that this reduces to 


det 


Ms 

NMi 


NMi 
Ms ’ 


(E25) 


which looks just like the three-mode instability equations with k ^ Nk. Therefore, we can read off the amplitude 
threshold immediately. Again, we see that the threshold is decreased by a factor of N compared to the three-mode 
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threshold. We expect the threshold energy to approximately scale as 


A 


2 

thr 


Nn J 4uj^k^ 


(E26) 


E.4. Decoupling of very different modes from collective sets 

In general, since all the mode parameters will be slightly different, our previous examples are a bit artificial. We 
now investigate the behavior when one mode begins to differ from the others. Consider the following characteristic 
equation, with N identical modes and one slightly different mode indicated by SM 




Ms Ml 
Ml Ms 


0 = det 

Mi Mi 
Mi Mi 

We can reduce this to 



0 = (det|Ms-M/|)^”^det 

Ms - Mi -SM 
NMi Ms + SM 

= (det \Ms — Mj\)^ ^ det M 5 -I- (A — 

l)Mi \ det 


Mi 

Mi 


Mi 

Mi 




(F27) 

Ms 

Mi 


Mi 

Ms+ SM 




(F28a) 

+ 

1 

{N - l)Mi)~'^ Mi + SM 

(F28b) 


As 6M —>■ 0, this reduces to a single collective set with N —>■ N + 1, as expected. We also note that this looks like the 
eigenvalues of a normal collective set with N modes and a new eigenvalue related to the different mode. Furthermore, if 
SM dominates the new eigenvalue, then we see that the different mode will “decouple” from the other modes. Clearly, 
there will be some threshold for how large SM needs to be before the different mode decouples, and that threshold 
will depend on the parent’s amplitude in a non-trivial way. We expect that a large parent amplitude Ag will support 
a larger SM before the mode decouples. 


F. SCALING OF PARAMETRIC INSTABILITY THRESHOLD AND WAVE BREAKING THRESHOLD 


Our calculations treat the system of modes as a set of global standing waves. ^__ 

parameter krfr 1, the wave will invert the stratification of the star and break ([Goodman & Dickson||1998 


& Ogilvie 


However, if a wave’s nonl i nearity 


Barker 


2010|. Because it does not reflect at turning points within the propagation cavity, such a wave is more 
ely treated as a traveling wave rather than a standing wave. Given that we specifically focus on parent 

< 1 ), we know that the parent is well described as a standing wave. 


appropriai 

waves below the wave breaking threshold {krf. 

Here we are intereste d in determining whether the same is true of the daughters, granddaughters, etc. 


As we describe in § 3.1 the parametric instability threshold scales as Athr oc u}°. This implies that each successive 
generation has a lower Athr and is therefore ever more susceptible to parametric instabilities. We show below that 
the energy above which a wave breaks also scales as Fbreak oc w®. Moreover, we find that Athr ^ Fbreak- This means 
that well before the daughters, granddaughters, etc. reach the wave breaking limit krfr ^ 1 , they will excite the next 
generation of modes through parametric instabilities. Although a mode is not necessarily limited to remain below its 
Ethr, we do not expect it to greatly exceed it either. This is because as a mode’s amplitude increases past its Fthr, its 
children grow at an ever faster rate and thereby limit how far their parent overshoots Ethr- While this issue requires 
further study, it suggests that our assumption that the modes are all global standing waves may be reasonable. 

We begin by calculating Ethr- For typical parameter values of a hot Jupiter system, Ethr is limited by the nonlinear 
detuning of the daughter modes rather than their linear damping (and similarly for granddaughters, etc.). To a 
first approximation, the detuning A is determined by half the frequency spacing between the daughter modes ujj2n. 
However, this assumes that the lowest Athr pairs are self-coupled modes. Because there is a distribution of mode 
frequencies slightly above and below half the parent frequency, there are always some mode pairs that happen to have 
A ^ uj/2n (Wu & Goldreich 2001). These are the pairs that minimize Ethr- Writing A = aujl2n, where a ^ 1 and 
using the expressions tor w, 7 , and k given in § | 2 . 2 [ we find that the threshold energy for self-coupled daughters is 


Ethr ^ 8 X 10 


-16 


Vo.oi/ 


P 

day 


-6 


An 


(FI) 


where a ^ 0.01 based on our three-mode network search results (cf. Figure!^. 

Now consider krfr- It is at its maximum near the inner turning point of the parent (where uj ~ N). This is because 
in the core of a solar model, k^ — AN/ujr is approxim ately constant and oc r~^ by flux conservatio n. Using the 
WKB relations given in Appendix A of WAQB (see also Goodman fc Dickson|199^ Ogilvie fc Lin|[2007 ), we find that 
the wave breaking condition max{A:j.G} = 1 for I = 2 modes corresponds to an energy 


Abrk — 3 X 10 


Uayy 
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Longer period modes break at lower amplitudes because they reach further into the core of the star. We thus see that 
both energies scale as w® and ifthr "C i?brk, as claimed. 


G. ESTIMATE OF THE LINEAR AND PARENT-DAUGHTER ORBITAL DECAY TIMESCALES 


The linear dissipation rate of individual resonant modes is — 2'yaE\i-n, where E\i-n is given by Equation (15). 


For the short periods that we consider, Aq, « ojal^Ua 7a- After summing over maw parents near the resonance, 
using the WKB estimates for the damping and forcing coefficients (Equations 8b andl^, and averaging according to 
Appendix [A} we find 


and 


~ 1.1 X 10 
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dayy 
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(G2) 


As we explain in § |4.2[ we can also estimate the nonlinear dissipation rate of networks consisting of only parents 
and daughters (but not granddaughters, etc.). This is because the dissipation in that case is dominated by the single 


daughter pair (/3,7) with the lowest instability threshold Ethr- As we show in Appendix 
a hot Jupiter system, the nonlinear equilibrium energy of such a daughter pair is Ep^^ ~ 


for the parameters of 
y/2Kap-y\EQ. The total 


dissipation rate of the system is approximately the dissipation due to these two daughters Ep.^ ~ 2 x 2'yp^^Ep^^. There 
is a small correction to this because the lowest Ethr daughters have slightly different parameters and therefore do not 
sit at exactly the same amplitudes. After accounting for this small correction and plugging in Equations [8bl and [TOl 
we find 


('’■£;)p-d — 2.0 X 10^ 
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and 


- 1.5 X 10® 
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Here we took = 1, which is representative of the typical lowest Ethr daughters for P > 2 days. We find good 
agreement between the parent-daughte r ne twork integrations that include many daughters and this analytic estimate 
(see circles and dashed curve in Figure [l0|). In the figure, we assume = 1 even for P <2 days. However, at these 
shorter periods, the available daughter modes are spaced further apart in frequency and the lowest Pthr pair may be 
pushed to > 1. This causes the small discrepancy between the circles and dashed curve at P < 2 days seen in 
Figure [TOl 


REFERENGES 

Aerts, C., Christensen-Dalsgaard, J., &; Kurtz, D. 2010, iSBN: 
978-1-4020-5178-4 

Alonso, R., Auvergne, M., Baglin, A., et al. 2008, A&A, 482, L21 
Bakos, G. A., Hartman, J. D., Torres, G., et al. 2012, AJ, 144, 19 
Barker, A. J. 2011, Monthly Notices of the Royal Astronomical 
Society, 414, 1365 

Barker, A. J., & Lithwick, Y. 2014, mnras, 437, 305 
Barker, A. J., Sz Ogilvie, G. 1. 2010, mnras, 404, 1849 
—. 2011, mnras, 417, 745 

Birkby, J. L., Cappetta, M., Gruz, P., et al. 2014, Monthly 
Notices of the Royal Astronomical Society, 440, 1470 
Bondarescu, R., Teukolsky, S. A., & Wasserman, 1. 2009, 

Phys. Rev. D, 79, 104003 

Brink, J., Teukolsky, S. A., Sz Wasserman, I. 2005, prd, 71, 064029 
Gandolfi, D., Parviainen, H., Deeg, H. J., et al. 2015, A&A, 576, 
All 

Gillon, M., Smalley, B., Hebb, L., et al. 2009, ASzA, 496, 259 
Gillon, M., Lanotte, A. A., Barman, T., et al. 2010, A&A, 511, A3 
Goldreich, P., Sz, Soter, S. 1966, Icarus, 5, 375 
Goodman, J., Sz Dickson, E. S. 1998, The Astrophysical Journal, 
507, 938 

Hebb, L., Collier-Cameron, A., Triaud, A. H. M. J., et al. 2010, 
ApJ, 708, 224 

Hellier, C., Anderson, D. R., Collier-Cameron, A., et al. 2011, 
ApJL, 730, L31 

Jackson, B., Barnes, R., Sz Greenberg, R. 2009, The Astrophysical 
Journal, 698, 1357 

Jackson, B., Greenberg, R., Sz Barnes, R. 2008, The Astrophysical 
Journal, 678, 1396 


Kumar, P., Sz Goodman, J. 1996, ApJ, 466, 946 
McQuillan, A., Mazeh, T., Sz Aigrain, S. 2013, ApJL, 775, Lll 
Meibom, S., Sz Mathieu, R. D. 2005, ApJ, 620, 970 
Mortier, A., Santos, N. C., Sousa, S. G., et al. 2013, ASzA, 558, 
A106 

Ogilvie, G. I. 2014, Annual Review of Astronomy and 
Astrophysics, 52, 171 

Ogilvie, G. I., Sc Lin, D. N. C. 2007, ApJ, 661, 1180 
O’Leary, R. M., Sz Burkart, J. 2014, mnras, 440, 3036 
Paxton, B. 2004, Publications of the Astronomical Society of the 
Pacific, 116, pp. 699 

Penev, K., Jackson, B., Spada, F., Sz Thom, N. 2012, The 
Astrophysical Journal, 751, 96 

Penev, K., Sz Sasselov, D. 2011, The Astrophysical Journal, 731, 
67 

Schenk, A. K., Arras, P., Flanagan, E. E., Teukolsky, S. A., Sz 
Wasserman, I. 2001, Phys. Rev. D, 65, 024001 
Smith, A. M. S., Anderson, D. R., Collier Cameron, A., et al. 
2012, AJ, 143, 81 

Storch, N. I., Sz Lai, D. 2014, mnras, 438, 1526 
Teitler, S., Sz Konigl, A. 2014, ApJ, 786, 139 

Terquem, C., Papaloizou, J. C. B., Nelson, R. P., Sz Lin, D. N. C. 
1998, ApJ, 502, 788 

Udry, S., Sz Santos, N. C. 2007, Annual Review of Astronomy and 
Astrophysics, 45, 397 

Valsecchi, F., Sz Rasio, F. A. 2014, ApJL, 787, L9 
Van Hoolst, T. 1994, aap, 286, 879 

Venumadhav, T., Zimmerman, A., & Hirata, C. M. 2014, The 
Astrophysical Journal, 781, 23 
Watson, C. A., Sz Marsh, T. R. 2010, MNRAS, 405, 2037 









ORBITAL DECAY OF HOT JUPITERS 


23 


Weinberg, N. N., Arras, R, Quataert, E., & Burkart, J. 2012, The Winn, J. N., & Fabrycky, D. C. 2015, Annual Review of 
Astrophysical Journal, 751, 136 Astronomy and Astrophysics, 53, 409 

Wu, Y., & Goldreich, P. 2001, ApJ, 546, 469 



